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Abstract 



We obtain analytic expressions for the third-order corrections due to the strong 
interaction Coulomb potential to the S'-wave Green function, energy levels and 
wave functions at the origin for arbitrary principal quantum number n. Together 
with the known non-Coulomb correction this results in the complete spectrum 
of ^-states up to order o^. The numerical impact of these corrections on the 
Upsilon spectrum and the top quark pair production cross section near threshold 
is estimated. 



1 Introduction 



Several years ago advances in non-relativistic effective theory made it possible to compute 
quarkonium properties at the next-to- next-to-leading order (NNLO). Since the pertur- 
bative approach assumes that the non-relativistic energy scale E ~ ma 2 s is larger than 
the strong interaction scale Aq CD , these computations apply to the lowest bottomonium 
states and heavy quark current spectral functions near threshold such as those relevant 
to top quark pair production. NNLO results have been obtained for matching coeffi- 
cients [HEllS], energy levels [H El El [7] , wave functions at the origin [El El [7] and spectral 
functions (6j [Tl [HJ, [9J, [TO] , principally for the S-wave states, which are the most important 
ones for applications [TT] . 

It was observed that the perturbative corrections are almost always very large. Al- 
though the origin of these large corrections can sometimes be understood as being due 
to mass renormalization [T2] or large logarithms, it is currently believed that a com- 
plete third (next-to-next-to-next-to-leading/NNNLO) order calculation is necessary to 
describe accurately even the case of top quark pair production near threshold. That this 
is a difficult undertaking is reflected by the fact that partial results at NNNLO exist for 
various quantities [131 H31 HSI HS1 H3 [^l HS1 12S1 12H 1221 12B], but only the S-wave ground 
state energy is currently fully known at third order [21], save for an unknown constant 
03 in the Coulomb potential. In this paper, we take one further step and compute the 
NNNLO corrections from the QCD Coulomb potential to the S'-wave energy levels and 
wave function at the origin for arbitrary principal quantum number n, and to the S'-wave 
Green function (spectral function) relevant to top quark pair production. Together with 
the known non-Coulomb energy level corrections [T8l [25] this determines the S-wave 
energy levels completely at third order for any n. The correction to the Green function 
and wave functions at the origin forms part of the complete third-order top quark pair 
production cross section. Another motivation for first concentrating on the Coulomb 
corrections is that the Schrodinger equation with the Coulomb potential can be solved 
numerically, and the result can be compared to the perturbative computation. This is 
no longer possible once the singular non-Coulomb potentials are included, in which case 
the perturbative approach is the only option. Comparing the perturbative approxima- 
tion to the numerical solution allows us to estimate the convergence of the successive 
appr oximat ions . 

2 Outline of the calculation 

The computation of Coulomb corrections can be phrased in the language of elementary 
quantum mechanics. We consider the Hamiltonian 
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where m denotes the heavy quark pole mass, and V the Coulomb potential of the strong 
force. With 
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the potential reads, up to the fourth order in the expansion in the strong coupling a 
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Here /3j are the coefficients of the QCD /^-function in the MS-scheme, defined with the 



convention da s /d\nfi 2 = —Y.Pn ( ^ +2 1 a f 
|26j . The constants ai, a 2 are given in 



= a s /(47r) such that /3 = — 4T F nf/3 

[3], the a^C^lnz/ 2 /q 2 term is from [13], while 
the three-loop constant term 03, currently unknown, is estimated to be 6240 (rif = 4) 
and 3840 (rif = 5) [27J. The group theory factors for SU(iV c ) gauge theory (N c = 3 in 
QCD) are C F = (N 2 - 1)/(2JV C ), C A = A^ c , Tp = 1/2, and n ; denotes the number of 
quarks whose masses are smaller than ma s and neglected (n/ = 4 for bottom systems, 
nf = 5 for top). a s refers to the strong coupling in the MS scheme at the renormalization 
scale /1, while v is a scale necessary to define the separation of potential and ultrasoft 
effects. The dependence on this factorization scale cancels in physical quantities when 
the ultrasoft correction is included. 

When 5V(q) = the spectrum of the Hamiltonian reproduces the well-known Bohr 
energy levels. This will be our zeroth order approximation. We consider quarkonium 
systems in which a s (ma s ) is small, and compute the spectrum in a perturbation expan- 
sion in a s . Hence, the coefficient of a™ +1 in 5V (q) is considered a perturbation of the 
nth order, and the third-order result we are aiming at requires up to three insertions of 
the a 2 potential, but only one insertion of the order potential. We shall focus on the 
5-wave Green function 



G(E) = (0\G(E)\0) = (0\ 
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where |0) denotes a position eigenstate with eigenvalue r = 0, and compute the matrix 
element of the right-hand side of (omitting the argument E of the Green function) 



G = Go — G SViG — GqSVzGo + G 5ViG 5ViGo 

- G 6V 3 G + 2G 5V 1 Go6V 2 G - GqSV^oSViGqSV^o. 
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Gq(E) denotes the zeroth-order Green function, and 8V n the nth order perturbation 
potential. The Green function G(E) has single poles at the exact S'-wave energy levels 



E = E n , such that 
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into this equation, we can determine 61,2,3 and by comparing the expansion of ([7j) 
in a s with © near E = E n . We recall that E$ = -m(a s C F ) 2 /(An 2 ), and |</>i 0) (0)| 2 = 
(masCp) 3 1 (87m 3 ) . The details of the calculation are too lengthy to be reproduced here. 
However, we sketch the computation of the threefold insertion (0\G 5ViGo5ViG 8ViG \0) 
of 5V\ in the Appendix. 

There is an equivalent quantum-field-theoretical description of the calculation, which 
is the appropriate one in the context of systematic higher-order computations of quarko- 
nium properties [11] . After all the quantum-mechanical description breaks down be- 
yond some accuracy, since (a) the system is sensitive to the short-distance quantum 
fluctuations, and (b) the restriction to the quark-antiquark sector of the Fock space is 
inadequate. Systematic calculations therefore start from QCD and obtain effective La- 
grangians by systematically integrating out short- distance fluctuations (Ax ~ 1/m) as 
well as all massless degrees of freedom down to the ultrasoft energy scale E ~ ma,, 
which characterizes the size of energy fluctuations in a quarkonium bound state. The 
result is an effective Lagrangian, in which the potentials appear as matching coefficients 
[28[ [29] . In particular the Coulomb potential is regarded as the matching coefficient of 
a relevant operator, and differs from the Wilson loop definition beginning at order oq. 
Consequently, the a A s lnct s term which appears in the static potential [30] is absent, and 
replaced by a logarithm of the factorization scale v [T3] as shown in (j3J). The computation 
of Coulomb corrections is based on the Lagrangian 
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where V(r) is the (colour-singlet) Coulomb potential, and where the non-relativistic two- 
component field ip(x) [x{%)) destroys (creates) a heavy quark (anti-quark). There are 
other terms in the effective Lagrangian that have to be included for a complete NNNLO 
calculation. The point we wish to make here is that these are well-understood, and hence 
there is a well-defined procedure to separate the complete calculation into several simpler 



3 



parts, of which the calculation of Coulomb corrections is one. An important feature of 
ffTUj) is that the leading-order Coulomb potential cannot be treated as a perturbation, 
and must be included in the zeroth order approximation together with the bilinear terms. 
(5 is essentially the propagator of this theory. We can also obtain G(E) in ([5]) directly 
by computing the two point function 

N C G(E) = i J d'xe^ (Q\T([^ X ](x)[x^}(0m, (11) 

with the Fock space vacuum state. This makes it clear that the Green function 
is closely related to inclusive heavy quark-anti-quark cross sections in e~e + collisions, 
where the case of the top quark is particularly interesting. 



3 Third-order corrections to bound state parame- 
ters 

In this section we present our result for the third-order correction from the Coulomb 
potential to the 5-wave energy levels and wave functions at the origin (squared) for 
arbitrary principal quantum number n. The corrections (notation as in ( 1H1I91) ) are pa- 
rameterized as 

e^ef + ef, (12) 

fi = ff + ff, (13) 

where C stands for the correction, when only the Coulomb potential is included in the 
effective Lagrangian. By definition l nC denotes all the remaining corrections, which 
originate from additional potentials and an ultrasoft non-potential interaction. Together 
with the non-Coulomb third-order correction e^ c to energy level [18, 25J, we obtain 
a complete result for the bound state energies to order mot 5 s for any n. For the wave 
function, however, f^ c is not yet known. We should emphasize that we define |?/> n (0)| 2 by 
the residue of the correlation function ( iTTj) of non-relativistic currents. The non-Coulomb 
corrections arise from potentials which cause short-distance singularities, such that /" 
are factorization scheme-dependent quantities. This scheme-dependence is canceled by 
short-distance coefficients, which we do not discuss here, but which are known to the 
same order as f™ c is known (i < 2). 

3.1 Energy levels 

The energy corrections from the Coulomb potential are organized as follows, 

e? = 4PoL + c E)1 , (14) 
4 = 12$ L 2 + L ( - 8$ + 4ft + 6Poc E>1 ) + c E>2 , (15) 
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where L = In {nn/ {mCpois)) and S a {n) = Sfc=iV^ a * s the harmonic sum. For later 
convenience we introduce the nested harmonic sums 
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In the following we omit the argument of harmonic sums which is always understood to 
be the principal quantum number n. The coefficients of the logarithmic terms are fixed 
by the renormalization group in terms of the /3-function and the ce,i from lower orders. 
The "non-trivial" information is encoded in the constants ce,%- The first and second 
order corrections, ce,x, c e,2 are known [U [5J, El E] 



c E ,i = 1a x + 4:SxP , 

ce,2 = a\ + 2a 2 + 4S x Pi + 4ai/3 



3S ± - 1 



2tt 2 



Sx 12Sx - 8 - - + 165 2 - 8nS 3 + — + 8<(3) 
nj 3 



(18) 



(19) 



where £(s) is the zeta-function, £(s) = Y^k=i k s - Our new result is the third-order 
correction to E n for arbitrary n, which reads 
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For completeness we also give the non-Coulomb correction. The first-order correction 
is only from the Coulomb potential, thus e™ c = 0. The second-order term was first 
obtained in [3], the third-order term for arbitrary n in [25]. The expressions are 
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where S 2 is the eigenvalue of the spin operator. For the spin-triplet (singlet) state 
S 2 = 2 (S 2 = 0). (The Coulomb potential is spin- independent, hence the ef do not 
depend on S 2 .) Furthermore L F (n) denotes the "Bethe logarithm", which must be 
evaluated numerically. For n — 1, 2, . . . we find 



L F (n) 



-81.5379, -37.671, -22.4818, -14.5326, -9.52642, -6.0222, 
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The first three numbers have been obtained in [14] . We have performed an independent 
calculation of the ultrasoft correction. 
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3.2 Wave functions at the origin 

The wave function corrections /j from the Coulomb potential are given by 
/f = 6#,L + £v,i, 
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The first and second order corrections are known [5j [7] 
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The first-order non-Coulomb correction vanishes, /" = 0, as for the energy-levels. 
Starting from second order ^(O)! 2 depends on the factorization scheme that separates 
the hard (relativistic) corrections from the low-energy corrections reproduced by the non- 
relativistic effective theory. The second-order term was first obtained in [5] for the spin- 
triplet states (5 2 = 2), but the result given there does not refer to the conventional MS 
subtraction scheme. The result of [5] including the hard correction has been reproduced 
in [7j, where the MS scheme was used for the individual contributions. Using these 
results (not printed in |7j) we find for the 5- wave spin-triplet non-Coulomb correction 
to the wave function at the origin squared in the MS subtraction scheme 
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where L(jj,f) = In (n/x f/(mCpa s )) and /// refers to the factorization scale. The \ij 
dependence cancels in the product C 2 |?/> n (0)| 2 , where C denotes the hard matching 
coefficient of the operator t/)'a % x [D 12] • (c l are the Pauli matrices.) The third-order 
coefficient f^ c is unknown. 
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4 Quarkonium masses 



In this section we compare the calculated energy levels with the bottomonium T(nS) 
masses. We also discuss the masses of the would-be toponium states, which are relevant 
to the top quark pair production cross section in high-energy e~e + collisions. 

First we give a numerical version of the general result for the 5-state energy levels 
for the spin-triplet case (5* 2 = 2), and rif = 4 (bottomonium) and rif = 5 (toponium). 
For the first three states n — 1,2,3 we obtain, for n 
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Here m b denotes the bottom quark pole mass, L = ln(n// / \m b C Va s (//))) , and we normal- 
ize the contribution from the unknown third-order constant in the Coulomb potential, 
a 3 , to the Pade estimate by defining a 3 = a 3 /a 3j p a ,i e . We have given the Coulomb and 
non- Coulomb corrections separately to emphasize the numerical dominance of the former 
(in the pole scheme) . The quarkonium masses are of course independent of the ultrasoft 
factorization scale u, but the separation into a Coulomb and non-Coulomb correction is 
not. The representations of the series above is given for v = m b Cpa s (fi)/n. We note 
that the Coulomb correction increases with n, while the non-Coulomb corrections be- 
come smaller. The reason for this is that the characteristic distance scale (the "Bohr 
radius") increases (r n ) ~ n, hence the relative effect of the short-range non-Coulomb in- 
teractions decreases for the excited states. The third-order result for n — 1 has already 
been obtained in [24], the other results are new. 



9 



For the spin-triplet toponium, n/ = 5, the series read 
4 



M t i(is) = 2m t - -m t a s 
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(36) 



where m t is top quark pole mass. 

The series coefficients are large and the series do not converge for bottomonium, 
presumably because the pole mass introduces a strong infrared renormalon divergence 
[3T], [32] , which is not present in the physical observable "quarkonium mass" itself [T2| 133] . 
It is therefore an advantage to use a better mass renormalization convention such as the 
potential-subtracted (PS) mass [12] . The relation to the pole mass appropriate to third- 
order calculations is given by 



"ips (/•*/) + 



/j,fC F a s 
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1 + ^(2(3 h + oi) + (a/3 2 l 2 + 2 (2oiA) + Px)h + a 2 ) 

(£)"( 8/5 ° 3/3 + 4 ( 3fll/3 o + ^A)i2 + 2(3a 2 /5 + 2a!/?! + &)Zi 
+ a 3 + 16tt 2 



(37) 



where h = ln(/i//i/) + 1, Z 2 = hi^/V/- 4 /) + 21n(/i//i/) + 2, Z 3 = ln 3 (/i//i/) + 3 ln 2 (/i//i/) + 
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61n(/i//i/) +6. Note that the PS mass is defined with the choice v = fif in the Coulomb 
potential V(q). The scale [if should be chosen of order ma s in order not to violate the 
power counting of the non-relativistic expansion, so the relation fl37|) is accurate to order 
ma b s just as the third-order bound state masses. Our standard choice is /if — 2 GeV for 
bottom and fif = 20 GeV for top. 



4.1 Bottomonium 

Before extracting quark masses and predicting bottomonium masses it is instructive 
to display the convergence of the expansions at the "natural" renormalization scale /z, 
where L = ln(n/i/ \M b CFCt s {^))) = 0. (Here M b refers to the bottom quark mass in the 
chosen scheme.) For Aqq^ = 290.4 MeV, and with 4-loop running of a s the "natural" 
scale is realized at /i = (2.02, 1.30, 1.03) GeV (for n = 1,2,3) with m b = 5 GeV, and 
H = (1.91, 1.23, 0.98) GeV with m 6iPS (2 GeV) = 4.6 GeV. For the numerical analysis we 
adopt a 3 = a^p a de = 6270 (S3 = 1). Eqs. (I31H33I) show that the precise value of 03 is not 
important as long as it is not very different from the Pade estimate. 

Using the relation fl37j) between the pole and PS mass we re-express Mrt n s) m terms 
of ttt^ps taking into account that fif/m b counts as one power of a s . With m b = 5 GeV, 
m^ps = rrib,ps(2 GeV) = 4.6 GeV we obtain 



M T(1S) = 2m b + e{ 0) 
= 2m b ^ s + e[°I s 



1 + 1.09 NLO + ( 1-42 + 0.36 nC ) N2LO + (2.29 + 0.28 nCj ^ ;j ( ( 



1 + 0.19 NLO + (0.07- 0.23 nC ) N2LO + (0.09- 0.19 nCyN ;j(( 



(38) 



M r{2S) = 2m b + Ef ] 



2^ife,ps + ^2jps 



M T(3S) = 2m b + Ei 0) 



2mf )> ps + E^p S 



1 + 1.91 NLO + (3.84 + 0.35 nC ) N2LO + (8.64 + 0.18 
1 + 0.26 NLO + (0.26 - 0.05 nC ) N2Lo + (0.37 - 0.03 n( 

1 + 2.69 NLO + (7.06 + 0.34 nC ) N2LQ + (20.10 - 0.03 nC 
1 + 0.25 NLO + (0.41 - 0.03„ C ) N2LO + (0.69 + O.OO.r^ 



N^LO 



N 3 LO 



(39) 



(40) 



with 

E n,PS = ^ + ( 41 ) 

This clearly shows that the series expansions are useless in the pole scheme, but the 
successive terms are (slowly) decreasing in the PS scheme for n — 1. For n > 1 we still 
observe that the transition to the PS scheme eliminates the huge correction from the 
Coulomb potential present in the pole scheme, yet the series coefficients are no longer 
converging. This is perhaps not surprising, because the scales are near or below 1 GeV 
for n > 1, and a perturbative treatment is simply no longer justified. 
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4.4 [ 

1 1.5 2 2.5 3 3.5 4 

H [GeV] 



Figure 1: The bottom PS mass, m^ps^ GeV), extracted from the experimental value 
Mxtis) = 9.460 GeV as a function of renormalization scale /x at LO (long dashes, black), 
NLO (long-short dashes, red), NNLO (short dashes, green) and NNNLO (solid, blue). 



We can therefore use the experimental T(1S) mass My( is ^ = 9.460 GeV to extract 
the bottom PS mass at NNNLO as was done in [31] at NNLO. In Figure [1] we show 
the extracted PS mass as a function of the renormalization scale \x at LO (long dashes, 
black), NLO (long-short dashes, red), NNLO (short dashes, green) and NNNLO (solid, 
blue). Varying \x from 1.25 GeV to 4 GeV (as done in [31]) the NNNLO correction is 
never larger than about 30 MeV and the error from the scale dependence is of similar size. 
We therefore assign a ±30 MeV error to m^pg from the truncation of the perturbative 
expansion. The uncertainty in a s [Mz) — 0.118 ± 0.003 results in a ±10 MeV error on 
m^ps- The largest uncertainty in the determination of the bottom quark mass from the 
T(1S) mass is then from non-perturbative effects. The perturbative approximation is 
justified when the ultrasoft scale mb(CV« s ) 2 ^> Aq C d, in which case the leading non- 
perturbative contributions is expressed in terms of the gluon condensate 

624vrm 6 (a s GG) 
425 (m^CVc^ 



SM r(is) = ,rX ru - (42) 



The numerical estimate is strongly dependent on the choice of scale in a s in the denomi- 
nator. Referring to |34J for a more detailed discussion of the non-perturbative correction, 
we obtain 

m fe ,ps(2 GeV) = (4.57 ± 0.03 pcrt . ± 0.01 Qs ± 0.07 non _ pcrt .) GeV (43) 

as the final result of this analysis. Because of the small third-order correction the bottom 
quark mass remains practically unchanged compared to the second-order analysis of 
[54"] . and so does the MS mass determined from (T43"]) . Further improvement of the mass 
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Figure 2: Predicted masses of the T(2S) and T(3S) as a function of the renormalization 
scale (jl. The lines refer to LO (long dashes, black), NLO (long-short dashes, red), 
NNLO (short dashes, green) and NNNLO (solid, blue). The widths of the bands for the 
experimental mass values ggerated. 



determination by this method requires a quantitative understanding of non-perturbative 
effects. 

Having determined m& ps, we are in the position to predict the masses of the excited S- 
level states at the third order. (An analysis of the complete spectrum at second order was 
performed in [37].) We focus on the spin-triplet states T(2S) and T(3S). The successive 
approximations up to the third order are shown in Figure [2] for m&^s = 4.57 GeV as 
a function of the renormalization scale. For /i between 2 GeV and 4 GeV it appears 
that the large third-order correction is welcome to bring the theoretical result closer to 
the observed masses. However, at the natural scales 1.23 GeV (n = 2) and 0.98 GeV 
(n = 3) the NNLO result agrees well with data and the NNNLO correction renders 
the prediction too large. As is apparent from the figure, the conclusion is that the 
perturbative computation of bottomonium masses is applicable only to the ground state, 
n = 1, while the excited states, involving larger distances, appear to be in the non- 
perturbative regime. It can be seen from (I39ll40p that the NNNLO term for n > 1 is 
dominated by the Coulomb correction. 

4.2 Toponium 

We briefly discuss the situation for the toponium masses. Toponium bound states do 
not exist in nature due to large decay width T t ~ 1.5 GeV of the top quark [38], however 
the remnant of the IS toponium state should be visible as an enhancement in the cross 
section e~e + — > tiX near threshold. The convergence of the series expansion for the 
toponium IS mass is therefore a good measure for the accuracy to which the top quark 
mass might be determined from this cross section [10] . 

The method suggested in [7J[29] relies on determining m t ,ps = ^t,ps(20 GeV) from the 
cross section measurement and obtaining the top quark MS mass from m ti ps, since both 
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relations are expected to be expressible in terms of well-behaved perturbative expansions. 
Adopting m 4j ps = 175 GeV, rif = 5, Aqq B = 208 MeV and the natural scale \i = 
32.6 GeV, where L = 0, we obtain 

Ma(is) = (350 + 0.85 + 0.05 - 0.13 + 0.01) GeV = 350.78 GeV. (44) 

The sum of the series varies only by about 60 MeV when the scale is taken between 15 
and 60 GeV, although the convergence is no longer satisfactory at the lower scale. The 
small uncertainty implies that m tt ps can be determined with little theoretical error from 
the cross section. The largest uncertainty in the determination of the top quark MS 
mass then results from the unknown four-loop coefficient in the relation between the MS 
mass and the pole mass, which is needed to convert m 4) ps to the MS mass, as already 
observed in [7]. This uncertainty is estimated to be around 100 MeV (see Table 2 of 
[TO]). Our conclusions are in good agreement with [39J, where one of us investigated 
the direct determination of the top quark MS mass from M t i(is) also using the NNNLO 
result for the IS energy level. We should emphasize that none of these estimates take into 
account electroweak corrections, which are non-negligible, and must be included in the 
mass relations and cross section prediction before a comparison with the experimental 
cross section can be attempted. 

5 Third-order Coulomb wave functions at the origin 
and Green function 

In this section we turn to the discussion of the 5*-wave Coulomb Green function, and 
to the wave function of the origin squared (residues of the Green function at the bound 
state poles). Since the third-order correction is not completely known, we include in 
this section only the Coulomb corrections as defined in Section [31 i.e. we also neglect the 
(known) non-Coulomb correction at second order. The series expansions seem to be out 
of control for the wave functions in the bottomonium system, hence we focus on the case 
of the top quark and set nf — 5. We also set = a 3t p a d e = 3840. 



5.1 Wave function at origin squared 



The numerical version of the general result for the 5-state Coulomb wave function at 
the origin squared reads, for n — 1, 2, 3, 



57T 



1 + a s ( - 0.4333 + 3.661 L) + a 2 s (5.832 — 5.112 L + 8.933 L' 



+ a£(- 13.73 + 6.446 In 



v 



m t C F a s 



+ 39.72 L - 22.91 L + 18.17Z/ 



(45) 



U fn\\ 2 { m tC F a s ) 3 



1 + a s ( - 0.1769 + 3.661 L) + c^(l0.19 - 3.861 L + 8.933 L 2 
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Figure 3: The Coulomb wave function at the origin squared for the ground state (n = 1) 
normalized by |^i°^(0)| 2 at //g = 32.6 GeV is shown as a function of the renormalization 
scale /i. The input parameters are m tPS (20GeV) = 175 GeV, v = m tj psCVa s (/i). The 
lines refer to LO (long dashes, black), NLO (long-short dashes, red), NNLO (short dashes, 
green) and NNNLO (solid, blue). 



+ 20.36 



6.446 In 



/ 2u 



,m t C F a t 



65.31 L — 19.09 L + 18.17i/ 



(46) 



1-03(0)1 



c 



(m t C F a s 

2167T 



1 + a s ( 0.07953 + 3.661 L ) + a] (13.27 - 2.609 L + 8.933 L 2 



+ (^(- 22.86 + 6.446 In ( 



/ 3u 



m t C F a, 



83.07 L - 15.27L 2 + 18.17 i/ 



(47) 



We show the successive approximations for n = 1 in Figure [31 As before, we re- 
expressed the expansion in terms of the PS mass, treating /i/ ~ m t a s . However, we note 
that contrary to the energy levels, the introduction of the PS mass does not change qual- 
itatively the behaviour of the expansion, neither would this be expected on theoretical 
grounds. Our reference top quark mass is m 4 PS (20GeV) = 175 GeV, and the ultrasoft 
factorization scale is chosen to be v = m 4i pgCjra s (/i) such that the corresponding loga- 
rithm in (I4"5l) vanishes. We also normalized the result to the LO wave function at the 
scale jiB = Tnt^sCFOts^B) = 32.6 GeV. It is clearly seen from the figure that the approx- 
imations converge, and that the inclusion of the new third-order correction stabilizes the 
prediction further, provided fi is larger than about 25 GeV. We find a similar behaviour 
for n = 2, 3, where, however, the enhancement of the wave function relative to leading 
order is about 50% [n = 2) and 100% (n = 3) rather than roughly 5% as at n = 1. 

It may be disconcerting that the perturbative expansion breaks down already at 
scales as large as 20 GeV, where the strong coupling is still small. A more detailed 
analysis shows that this early breakdown is caused primarily by the a s (a s /3o In q 2 ) n terms 
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in the Coulomb potential. We also note that at /i = /jb, where L = 0, the series 
expansions shown above are sign-alternating in contrast to the corresponding expressions 
for the energy levels, which exhibit fixed-sign behaviour. It is a general fact that for 
sign-alternating series of a certain regularity, the convergence of the expansion is much 
improved by choosing a larger scale |40j . since this renders the series coefficients and the 
expansion parameter a s smaller. This explains the stability seen in the figure towards 
scales larger than the natural scale fj, = hb, and suggests that an error estimate for 
^(O)! 2 from varying the scale between /is/2 and 2fi B m ay be misleadingly large. We 
shall see in the following subsection that this is indeed the case. 

5.2 Green function 

In addition to the energy levels and wave functions at the origin we have also computed 
the full .S-wave Green function up to the third order in the presence of the Coulomb 
potential as described in Section [2j The result is expressed in terms of multiple sums 
that can be evaluated numerically only. 

The Coulomb Green function plays an important role in the calculation of inclusive 
top quark pair production e~e + — > tiX near threshold [UJ, since the bulk of the cross 
section is given by 



where e t = 2/3 is the top quark electric charge, T t the top quark width, and az accounts 
for the vector coupling to the Z boson. The convergence of the perturbative approxi- 
mation up to the second order has been the subject of many investigations several years 
ago (see the review [TU]). 

We are now in the position to extend this investigation to the third order as far as the 
corrections from the Coulomb potential are concerned. We computed the perturbative 
expansion of the Green function given in ([6]). This strict expansion is never a good 
approximation, because it contains terms of the form 



which originate from the expansion of ([7]) around E^' rather than the true pole position. 
These terms become numerically large near EttE£>\ but they can be summed by adding 
the exact pole structure and subtracting the expanded structure to the appropriate order, 
see [7J for the corresponding expressions. In the following, "perturbative approximation" 
means that this resummation is included. Alternatively, we computed the Green function 
([5]) numerically by solving the Schrodinger equation with the Coulomb potential (j3J) ex- 
actly, following the method described in [32] • We shall refer to this as the "exact result". 
The exact result contains an arbitrary number of insertions of the perturbation potentials 




(48) 



(49) 
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8Vi, SV2, SV3. The leading difference to the third-order perturbative approximation con- 
sists of fourth-order terms of the form G 5ViGoSViGo5ViG 5ViG , Gq5ViGq8ViG 5V2G 
etc. Comparing the two results we obtain an estimate of the importance of these multiple 
insertions and of the convergence of the perturbative approximation. 

For the following numerical study we assume A^" 5) = 0.208 GeV (a s (M z ) = 0.118), 
four-loop evolution of a s , T t = 1.5GeV, and v = 20 GeV. It is crucial for an accurate 
prediction of the cross section not to use the top quark pole mass as an input parameter 
[3 Hn]. Our result is presented in terms of the top-quark potential- subtracted mass 
^t,ps(20 GeV) = 175 GeV defined in (|37j) . The conversion to the top quark MS mass 
is given in [10J. We implement the PS scheme by working with an order- dependent 
pole mass. That is, in the N fc LO perturbative approximation we compute the top quark 
pole mass rnf^ from m t ,ps(20GeV) = 175 GeV according to fl37|) including the terms 
up to order yu^aj +1 , and use m[ as an input to the Green function. For instance, for 
renormalization scale fi = 30 GeV we obtain mf = (176.21, 176.56, 176.74, 176.87) GeV 
for k = 0, 1, 2, 3. The energy argument of the Green function is then 

E = ^~s- 2m ( j k) (50) 

with \fs the center-of-mass energy of the e~e + collision. 

The upper plot in Figure H] shows the convergence of the successive perturbative 
approximations to the exact result for the Green function (cross section) at the renor- 
malization scale /i = 30 GeV. The location of the "peak" position is indeed stable under 
the inclusion of higher-order corrections as expected in the PS scheme. On the other 
hand, the corrections to the magnitude of the cross section near the peak are significant, 
decreasing from about 20% at NLO to about 7% at NNNLO. The corrections alternate 
in sign as expected from the behaviour of the series expansion of |t/>i(0)| 2 . An important 
observation is that the third-order perturbative approximation coincides with the exact 
result within 1%. Hence the higher-order insertions of the perturbation potentials are 
negligible. The convergence of the approximations then suggests that the residual error 
from yet higher perturbation Coulomb potentials SV^ etc. should be less than 5%. 

In the lower plot of Figure H] we display the renormalization scale dependence of the 
third-order perturbative approximation (solid lines, \i = (60, 30, 15) GeV from top to 
bottom). It is immediately apparent that the scale dependence is very small from 30 
to 60 GeV, but the result for \i — 15 GeV is far away. We can trace this anomalous 
behaviour directly to the breakdown of the perturbative expansion for |?/>i(0)| 2 at scales 
below 25 GeV discussed in the previous subsection and displayed in Figure [31 The exact 
result (dashed line) does not exhibit this behaviour for /1 = 15 GeV, hence we conclude 
that the multiple insertions of the perturbation potentials become large at small scales 
and destroy the agreement of the perturbative and exact result. Indeed, we find that the 
series of multiple insertions is very slowly converging at small scales. We therefore learn 
the important lesson that the "correct" choice of scale in the perturbative approach is 
fi > 25 GeV, while choosing smaller scales would lead to misleadingly large uncertainties. 
The lower plot of Figure 0] indicates that the scale dependence is less than 5%, similar 
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Figure 4: Top quark pair production cross section (Coulomb corrections only) for m 4i ps = 
175 GeV, r t = 1.5 GeV. Upper panel: successive approximations up to the third order 
for n = 30 GeV. Lower panel: Scale dependence of the third-order approximation. See 
text for further explanation. 
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in size to the truncation error estimated above. This discussion does not include the 
unknown third-order non-Coulomb corrections, but it lends support to the optimistic 
interpretation that the magnitude of the ti threshold cross section can eventually be 
computed with an accuracy of a few percent. 

6 Conclusion 

We computed the third-order corrections from the strong interaction Coulomb potential 
to the energy levels and the wave functions at the origin of the S-waves bound states, and 
to the expectation value of the S'-wave Green function operator in the \r = 0) state. We 
view this as a first step towards a complete calculation of the third-order top quark pair 
production cross section e~e + — > tiX near threshold. It also completes the expression 
for the S'-wave quarkonium masses with accuracy ma\ for arbitrary principal quantum 
number n, except for the unknown constant as in the Coulomb potential. 

We updated the determination of the bottom quark mass from the mass of the T(1S) 
and obtain 

m fe ,ps(2 GeV) = (4.57 ± 0.03 pcrt . ± 0.01 Qs ± 0.07 non _ pert .) GeV (51) 

for the PS mass, with almost no modification compared to the second-order analysis 
|34j . Our numerical study of the Coulomb corrections to the top quark pair produc- 
tion cross section near threshold shows that the perturbative approach (mandatory once 
non-Coulomb corrections are included) works and led to the conclusion that the resid- 
ual uncertainty from the Coulomb corrections is less than 5%. This lends support to 
the optimistic interpretation that the magnitude of the ti threshold cross section can 
eventually be computed with an accuracy of a few percent. 

Note added 

During the preparation of this paper Penin, Smirnov and Steinhauser [25] have obtained 
results for the Coulomb correction to the second and third energy level and wave function 
at the origin (n = 2, 3), which agree with our result for general n. We thank the authors 
for communicating and comparing their results prior to publication. 
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Appendix 

The most difficult part of the third-order calculation of the Green function is the threefold 
insertion (0\G 8ViG 5ViG 5VxG \0) of the first-order perturbation potential 



4:7iC F a s a s \ fi 21 
a x + fj In — 
q 2 



(52) 



q 2 Air 

The matrix element is ultraviolet and infrared finite, so the computation can be done in 
three dimensions. Since the external state is a position eigenstate, we Fourier transform 
to position space, where the potential has a 1/r and a ln(r)/r term. We can generate 
these terms by working with 

= — , ^-j-, (rT^ (53) 



47rr(l + 2Ui) COs(7TMj) 

and by taking the zeroth and first derivative with respect to the at Ui = 0. The 
threefold insertion of the generating potential reads 

3 

J 3 = / nrf 3 nGo(0,r 1 ;E)^ 1 (r 1 )Go(r 1 ,r 2 ;E)^ 2 (r 2 )Go(r 2 ,r3;^)^3(r3)Go(r3,0;E) 

, r 2ui+l r 2« 2 +l r 2w 3 +l 

= / dr±dr 2 dr 3 — ^— -— 2 — - — ^— - 

J T(l + 2u\) cos(ttui) T(l + 2u 2 ) cos(7ru 2 ) T(l + 2m 3 ) cos(7to 3 ) 

x{g o (0, n, E)Gt 0) (ri, r 2 , £)G# =0) (r 2 , r 3 , £)G (0, r 3 , £)}. (54) 
For the zeroth-order S'-wave Coulomb Green functions we use the representations [43j 144] 

G? (0,ri,S) = —^ etmvn I dte 2mvrit (^-j^J , (55) 



2tt (n+ l)(n+ 1 - A) 



(56) 



with u = y (E + ze)/m and A = %C F a s /{2v). The L^(x) are the Laguerre polynomials. 
The integrals over can now be factorized into two functions H(u,n) and K(u,n,j), 
and we obtain 

t _ f m V V V H( Ul ,n)K(u 2 ,n,j)H(u 3 ,j) 
3 UJ ^^( n + l)(n+l-A)0- + !)(; + ! -A)' 1 ] 

where (defining s = —2imvr) 

- ^ (iiv)"r*( i ^) A r' i -- <i+,,,s2 " +i ^ , w- < 58 > 

K(u,n,j) = 



r(l + 2m) cos(7to) Am 2 v 2 \Am 2 v 2 

20 



J^dss 2u+1 e- s L^\s)Lf\s). (59) 



Performing the integrations in H(u, ft) (substituting x — 1/(1 + 1) in the t-integral), 
we find 



H(u,n) = (-AmE) 



_Jl + 2«)(n+l) 



COS (TTU) 



dtr x (i + t) 



\-2-2u 



2 F 1 -n, 2 + 2m; 2; 



1 



1 + t 



(n + l)r(l-A) " (-l) k n\ T{2 + k + 2u)T{l + k + 2u) 
~~ cos(7ru)(-AmE) u ^ Q k\(n-k)\T(l + 2u)T(k + 2)T(2 + k + 2u- A)' ^ ^ 

The sum can be expressed in terms of the hypergeometric function sF2( — n, 2 + 2m, 1 + 
2m; 2, 2 + 2m — A; 1), but it is simpler to perform the expansion in the generating variable 
u directly. We need the first two terms in the expansion. With the help of the generating 
sum 

'-l) k n\ T(l + k + a) r(l + a)F(l - a + n - A) 



E 



we find 



H(0,n] 
H'(0,n 
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(61) 



(62) 
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k J k\(n-k)\T(2 + k-X) 

+ 2^(1 + k) + 2^(2 + fc) - 2^(2 + fc - A) J 
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27s - In(-AmE) 
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- \n(-AmE) -2^ E + *(n + 1 - A) 
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*(1 - A) -*(n + 2- A) 



(63) 



(ft + 1) 

where ^{z) denotes Euler's Psi-function, and 7# = 0.577216 ... is Euler's constant. This 
result has already been obtained in [7J. 

Similarly, for K(u,n,j) we need the expansion up to the first order in u. Here we 
obtain 

ft + 1 



K(0,n,j] 



1 



Am 2 v 2 



I+{n+l)6 nj {2 lE -ln{-imE)) 



with 



ds2s\n(s)e- s L^\s)Lf\s). 



(64) 
(65) 

(66) 



To solve the integral /, the Laguerre polynomials are expressed through their generating 
functions. Then with 



_ z u 

e !-« 



(67) 
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and 



ds 2s ln(s)e 



s v _ s w q 

X—v p 1 — w £ 
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the integral J is expressed as 
1 1 <9 n d j 



[VW 



n\ j\ dv n dwi 



2 -1+7^ + ln 



l—vw 
(v-l)(w-l) 



2 + 2(l + n)^(l + n) 

2 min(n, j) + 1 
|j - n| 

Hence the final result for K'(0,n,j) reads 



(t>w — l) 2 

if n = j 
if n / j 



t)=u)=0 



K'(0,n,j) 



2 + (n + 1){2[7 jB + ^(1 + n)] - ln(-4m£?)} if 



n = j 



Am?v 2 



_ 2 min(n, j) + 1 
\j ~n\ 



if n ^ j 



(68) 



(69) 



(70) 



At this point, we have expressed the threefold insertion of the perturbation potential 
in terms of a doubly infinite sum involving Euler's Psi-function, see ( |57l) . These sums 
converge rapidly when the energy argument is evaluated along a line parallel to the real 
required in the calculation of the top quark pair production cross section. In 
order to obtain the third-order correction to the energy levels and wave functions at 
the origin in Section [31 we extract analytically the pole part of J3 when E approaches 
the leading-order S"-wave bound state energies -E^ - 1 , which correspond to A = n. This 
results in multiple sums, which, after a tedious reduction, can all be expressed in terms 
of zeta-f unctions and nested harmonic sums. 
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